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In  performing  engineering  or  scientific  data  analysis  or  computa¬ 
tions  it  frequently  becomes  necessary  to  examine  data  which  is  a  single 
valued  function  of  two  independent  variables.  One  convenient  method  of 
displaying  this  type  of  data  is  with  contour  plots.  IMs  report  describes 
an  efficient  algorithm  for  construction  of  contour  lines  and  the  imple¬ 
mentation  of  this  algorithm  as  a  FORTRAN  IV  subroutine,  CONTUR. 
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I.  INTRODUCTION 


In  performing  engineering  or  scientific  data  analysis  it  frequently 
becomes  necessary  to  examine  data  which  is  a  single  valued  function 
of  two  independent  variables.  A  common  and  useful  technique  for  dis¬ 
playing  such  data  is  through  the  use  of  contour  plots.  When  three 
independent  variables  'are  involved  any  method  of  graphic  display  is 
cumbersome,  but  the  plotting  of  contours  in  two  dimensions  for  several 
values  of  the  third  independent  variable  may  be  the  most  practical 
alternative. 

Organizations  heavily  involved  in  scientific  computation  utilizing 
digital  computers  frequently  need  to  reduce  results  into  an  easily 
comprehendable  form  such  as  contour  plots.  Accordingly,  it  Is  highly 
desirable  that  an  easy  to  use  subroutine  for  contour  plotting  be  available 
for  the  organization's  computer  users.  C0N1UR  is  such  a  subroutine, 
written  in  FORTRAN  IV,  and  hence,  is  compatible  with  many  digital  com¬ 
puters  in  use  today.  The  subroutine  described  herein  was  designed  to 
work  In  conjunction  with  the  California  Computer  Products,  model  780 
digital,  incremental  plotting  system  and  the  associated  plotting  sub¬ 
routines  in  use  at  BRL,  However,  with  simple  modifications,  CONTUR  may 
be  used  with  other  forms  of  graphic  display  equipment. 

II .  ALGORITHM 

The  data  for  which  contours  are  to  be  drawn  is  assumed  to  be  a 
discrete  tabulation  of  the  single  valued  function 

Z  =  f(x,y)  (1) 

for  x,y,  in  the  range  over  which  contours  are  desired.  For  a  fixed  Z, 
Z=2o>  Eq,  (l)  may  be  written 

Y  =  g(x,z0).  (2) 
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In  this  form  the  curve  is  called  a  contour  and  in  general  a  different 
contour  would  occur  for  each  value  of  Z  .  Usually  the  function,  f(x,y) 
is  not  known,  the  data  arising  either  from  experiment  or  by  numerical 
approximation  techniques.  Hence,  the  explicit  expression  as  a  func¬ 
tion  of  x  and  Z0,  Eq.  is  not  available  and  a  numerical  procedure 
for  determining  the  contours  is  necessary. 

The  algorithm  described  below  represents  a  significant  simplifica¬ 
tion  of  the  algorithm  described  by  James  Downing  [l ] . 

The  algorithm  is  derived  by  focusing  attention  on  four  adjacent 
data  points  Z.  Z  ,  Z.  .  and  Z  .  .  where  the  coresponding 
independent  variables  have  the  values  (X^  Y^),  Y^)  etc.  Assuming 

the  data  contains  I  points  in  the  x  direction  and  J  points  in  the  Y 
direction,  the  algorithm  must  be  applied  to  N  =  (I-l)(J-l)  cells. 

Within  such  a  cell,  Figure  1,  the  center  point  is  located  and  assigned 
a  Z  value  equal  to  the  average  of  the  four  Z  values  at  the  corners.  These 
five  points  are  then  connected  with  line  segments  which  are  in  turn 
numbered  one  through  eight  in  a  clockwise  direction. 


Figure  1 .  A  Typical  Cell . 
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Each  segment  is  then  tested  to  see  if  the  required  contour  inter¬ 
sects  with  it,  in  the  following  manner.  Starting  with  segment  one,  the 
contour  value  Z0  is  subtracted  from  the  end  points. 


T  =  Z  -  Z 
2  i+l,J  o 


(3) 


If  the  quantity 


(M 


is  greater  than  zero,  the  entire  segment  is  either  above  or  below  7  , 

if  A  equals  zero,  either  7  or  Z  is  equal  to  7.  ,  and  if  A  is 

1»J  J  0 

less  than  zero  the  contour  intersects  the  segment.  In  this  last  case 
the  point  of  intersection,  xq  is  found  by  linear  interpolation  ( se»* 
Figure  2)  with  xQ  given  by 


=  (Z 


Zi,J)Ui+l 


'  V/(Vl,j 


i.j) 


The  x  and  y  values  of  this  intersection  are  then  stored  ir  temporary 
arrays  PX  and  PY. 


11 


The  procedure  > g  then  repeated  for  segments  two  through  eight.  When  seg¬ 
ment  eight  is  completed,  the  points  stored  in  PX  and  FY  are  plotted  and 
the  next  set  of  points  are  considered. 

Before  the  ordered  pairs  (PX,  FY)  can  be  plotted  successfully  there 
are  several  condition^  which  must  be  tested  for  and  if  present,  properly 
handled,  (l)  If  all  four  of  the  cell's  corner  points  are  equal  to  7  ,  no 
points  should  be  plotted.  (2)  When  the  contour  intersects  segment 
eight,  the  PX  and  FY  arrays  must  be  reordered.  The  reason  for  this 
becomes  obvious  when  one  remembers  that  the  segments  are  tested  in  a 
clockwise  direction.  For  instance,  assume  C0N1UR  finds  intersections 
on  segments  one,  seven  and  eight.  Plotting  these  points  as  originally 
stored  would  result  in  an  extraneous  line  being  drawn.  See  Figure  3* 

By  simply  rearranging  the  points  so  that  they  are  stored  Beven,  eight, 
one,  the  correct  contour  is  drawn. 


Extraneous  Line 
Correct  Line 


(3)  Provision  is  also  made  for  the  case  where  two  contours  of  the  same 
value  pass  through  the  cell.  This  occurs  only  when  two  opposite  Z 
values  are  greater  than  and  the  other  two  points  are  less  than  Zq. 

By  noting  if  the  center  point,  Z2,  is  greater  than  or  less  than  ZQ,  the 
paths  taken  by  the  contours  are  specifically  known  and  are  plotted  as  a 
special  case.  See  Figure  k. 
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Ill .  THE  SUBROUTINE 


CONTUR  is  accessed  through  the  statement 

CALL  CONIUR  (Z,  X,Y,  IS,  IY,  DZ,NZ,  IZ) . 


Z,X,Y  are  the  arrays  containing  the  values  7  ,  X  and  Y  ,  respectively. 

IX  and  IY  are  the  number  of  points  in  the  X  dimension  and  Y  dimension. 
The  subroutine  requires  that  Z  be  of  dimensions  (IZ,IY).  DZ^NZ)  is  a 
one  dimensional  array  containing  the  7  values  at  which  contours  are 
desired.  NZ  is  the  number  of  these  values.  The  declared  number  of 
rows  in  the  Z  array  is  IZ. 


Since  the  subroutine  uses  just  four  data  points  at  a  ti  e,  requiring 
no  knowledge  of  where  it  has  been  or  ..ere  it  is  going, enormous  amounts 
of  data  can  be  handled  by  reloading  the  Z  array  and  calling  CONTUR 
several  times  with  different  portions  of  the  data. 
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The  computer  time  required  by  CON1UR  depends  on  the  size  of  the  7. 
array,  the  number  of  values  at  which  contours  are  desired  and  the 
smoothness  or  irregularity  of  the  data,  with  time  increasing  for  large 
arrays,  large  numbers  of  contour  values  and  irregular  data.  For  some 
typical  times  see  the  examples  of  contour  output. 

i 

In  order  to  keep  the  subroutine  as  efficient  and  machine  independent 
as  possible,  no  labeling  of  contours  is  done,  nor  are  any  borders  or 
titles  plotted  in  CONTUR.  The  user  must  initialize  the  plot  routines 
and  set  scales  prior  to  calling  CONTUR.  PLTCCD  is  a  predefined  sub¬ 
routine  on  the  Ballistic  Research  Laboratories  BRLESC  computers  that 
generates  input  data  for  the  Cal  Comp  780  digital,  incremental  plotting 
system,  and  must  be  replaced  for  use  of  CONUJR  on  other  computer  systems. 


BRLESC  users  should  note  that  the  positioning  on  the  plotter  page 
and  scales  used  by  C0N1UR  are  determined  by  the  latest  reference  to  PLTCCS. 
[l]  Thus,  it  may  be  necessary  to  reset  the  plotting  scales  before  cal¬ 


ling  CONTUR. 


IV.  CONTUR  INRJT  VARIABLES 


Z(IZ,IY) 

X(IX) 

Y(IY) 

IX 

IY 

DZ(N2) 

NZ 

n 


is  a  two  dimensional  a-ray  containing  the  functional 
values  of  the  data. 

is  a  one  dimensional  array  containing  the  values 
of  one  of  the  independent  variables. 

is  a  one  dimensional  array  containing  the  values 
of  the  other  independent  variable. 

is  the  number  of  elements  in  the  X  array. 

is  the  number  of  elements  in  the  Y  array. 

is  a  one  dimensional  array  containing  the  7  values 
at  which  contours  are  desired. 

is  the  number  cf  elements  in  the  IC  array. 

the  number  of  declared  rows  in  the  7  array. 


V.  FLTCCD  INKJT  VARIABLES 


N=l,  L=0 

px(i),pr(j). 

K 

M=0 

M=1 


signifies  that  a  line  plot  Is  to  be  drawn. 

are  arrays  containing  the  X  and  Y  coordinates  to  be 
plotted.  The  First  point  plotted  is  Px(l),PY(J). 

Is  'the  number  of  data  pairs  to  be  plotted 

causes  the  subroutine  to  start  a  new  curve  vlth  the 
point  (PX(I),  FY(J)). 

causes  the  subroutine  to  continue  the  curve  plotted 
by  the  previous  FLTCCD  entry. 
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APPENDIX  A 


Examples  of  CONTUR  Output 


Figure  A-i  is  a  three  dimensional  graph  of  the  function 

Z  -  SIN  (x  +  y)/( 1  +  (x  -  y)2)  A-l 

as  plotted  by  the  subroutine  GRAF3D^3\  Although  this  plot  is  interesting 
and  demonstrates  general  trends,  it  is  virtually  impossible  to  retrieve 
any  useful  quantitative  information  from  it.  Figure  A-2  is  a  contour 
plot  as  drawn  by  CONTUR  of  the  same  data.  The  contour  lines  are  at 
values  of  .1,  .4,  .6,  and  .9.  The  E  array  contains  36OO  points  and 
CONTUR  required  14.4  seconds  on  the  BRLESC  II  computer  to  generate  the 
curves . 

Figures  A-3  and  A-4  represent  experimental  data.  Again  the  contour 
plot  is  the  more  analytically  useful,  although  not  as  esthetically 
pleasing  as  the  3-D  plot.  The  PR01  subroutine  (4)  was  used  in  this 
case  to  generate  the  three  dimensional  plot. 

Figures  A- 5, 6, 7  are  included  to  demonstrate  the  results  of  allowing 
the  data  gride  to  become  too  @oarse.  Figure  A-5  is  the  3-D  representation 
of 

Z  =  |SIN  ^  x2  +  2)/(Vx2  +  y2  I  -20  <  x,  y  <  20. 

«/ 

Both  Figures  A -6  and  A-7  are  contour  plots  of  the  above  function  with 
contour  lines  drawn  at  Z  values  of  .1,  .4,  .6,  and  .9.  In  Figure  A-6 
the  grid  contains  10000  points  and  the  representation  is  accurate.  The 
grid  in  Figure  A-7  contains  only  2500  points  and  the  interpolation 
scheme  is  no  longer  sufficiently  accurate  to  portray  the  function  cor¬ 
rectly.  The  BRLESC  I  computer  required  60.6  seconds  to  generate  A-6 
and  19.2  seconds  for  A-7. 


IT 


A -2.  Contours  of  Z  =  SIN  (x  ;  ;/(l  +  (x  -  Y)‘ ) 
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A-3.  3-D  Plot  of  Experimental  Data 
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APPENDIX  B 


30,33 


YES 


Find  X  &  Y  values  at  intersection 
interpolating  if  necessary.  Store 
in  PX  &  FY  and  set  flag. 


39 


Does  the  contour  intersect  line 
segment  4? 


40 


YES 


Find  X  &  Y  values  at  intersection 
interpolating  if  necessary.  Store 
in  PX  &  FY  and  set  flag. 
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Does  the  contour  intersect 
line  segment  5? 

NO 

28 
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SUBROUTINE  CCNTUK  (Z,  X,  Y,  IX,  IY,  OZ,  NZ ,  IZ) 

DIMEN  SIC.)  X(  I  X  )  •  Y  (  IY),Z(  I  Z  •  I  Y )  »PX(8)  »  P  Y  I  8  )  ,KCHK ( 8 ) , 07  I NZ  ) 

THIS  SUBROUTINE  PLOTS  NZ  CONTOURS  AT  OZ  VALUES. 

X  AND  Y  ARE  ONE  DIMENSIONAL  ARRAYS  wF  LENGTH  IX  ANO  IY,  RESPECTIVELY 
Z  IS  A  TWO  UIMFNSI CNAL  ARRAY  OF  SIZE  (IX, IY). 

OZ  IS  A  ONE  0 1  MENS  I C  'AL  ARRAY  .F  LENGTH  NZ  IN  WHICH  THF  Z  VALUES  AT 
WHICH  CONTOURS  ARF  DESIRED  ARC  PLACED. 

THIS  VERSI  IN  OF  CONTUR  WAS  COMPLETED  IN  FERRUARY  1973. 

ic=: 

NNX  = I  X-  1 

NNY  =  IY-1 

DO  2  J=  i  »NNY 

DO  1  1=1,  NNX 

DO  100  ICON T  =  1 , NZ 

ZG  =  DZ (  IC  .NT  ) 

IF  ALL  FOUR,  DATA  POINTS  ARE  EQUAL  TO  ZO,  DO  NOT  PLOT  ANY  LINES  FOR 
THIS  CELL. 

IF(Z(  I  ,  J  )  .EU.ZC.AND.ZI  I«-l,.i).EO.ZO.AND.Z!  I  , J  +  l )  .  EQ.  ZO  •  ANO. 

1Z(  I  +  1,J*1).EC.ZC)  GOTO  ICC 

TEST  SEGMENT  1  FOR  AN  INTERSECTION  WITH  THF  CONTOUR  LINF. 

T  1  =  Z  I  I, J  )-Z0 

T2=Z< I+1,J)-ZC 

D=T1*T2 

IF(  0) 10,  11,  19 

10  IC=IC+1 

PX(ICI=-T1*IX(I+1)-X(I))/(Z(I+1,J)-ZII,J))+X(I) 

PY(IC)=Y(J) 

KCHK ( 1)=1 
GOTO  19 

11  IF(Tl.NF.O.C)  GOT:  13 
IC=ICM 

PXI  IC)  =  X( I ) 

PY (  IC  )  =Y ( J  ) 

KCHK (  1 »  =  1 
I F ( T  2 ) 19, 13,  19 
13  IC=IC+1 

PXI  IC)  XI  1  +  1  ) 

PY(  IC  )  =Y ( J ) 

KCHK (1)^1 

TEST  SFGMLNT  2  FOR  AN  INTERSECTION  WIJH  THE  CONTOUR  LINE. 

19  T3=.25*(Z(I,J)+Z(I+l,J)+Z(l,J+l)+Z(I+l»J+l> )-ZOt 
D=T 1*  T  3 

x 2=  ( x i  i  +  n+xm  >*.5 
Y  2=  (  Y  (  J  ♦  1 )  ♦  Y  (  J  ) )*.5 
Z2=T3+Z 
I  F ( D ) 20,  23,  29 

20  IC=lC+l 
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px(  ic ) *- t i* c  x2-xn )  )/iZ2-zci . jn*xm 

PY(  IC)=-TI*(  Y2-YU)  l/li2-Zll,J)  )  ♦  Y  C  J) 

KCHK ( 2 ) = 1 
GCTC  29 

23  IF(T3.Nc.0.C)  GCTC  29 
IC=IC*1 
PX(  I  C  I  =  X  2 
P  Y  (  I C  )  =  Y  2 
KCHK I  2  )  =  1 

TEST  SEGUE  IT  3  FOR  AM  INTERSECTION  WITH  THE  CONTOUR  LINE. 

29  T2  =  Z (  i« J+1I-ZO 
D=  T 1  *  T  2 

I  FI  0) 30, 3  3, 3  9 

30  IC=IC*1 

PX (  IC  )  =  X (  I  ) 

PYI  IC ) =- T 1* I  Y  I  J  ♦  1  l-Y(J)  )  /  (ZII  » J* 1 ) -Z ( I  ,  J )  )  ♦  Y  I  J  ) 

KCHK! 3 )  =  l 
GOTO  39 

33  I  r-  <  T2.NE.C.C  )  GOTO  39 
I  C=  I C ♦ 1 
PXI  T C  )  -=  X  (  I  ) 

PYI  IC  )=Y ( J  +  I  I 
KCHK l 3  I  =  I 

TEST  GFGMENT  9  FOR  AM  INTERSECTION  WITH  THE  CONTOUR  LIME. 

39  T  1  =  Z  l  I,  J  +  D-ZO 
U=T  1*T3 

IFIO) A 0,49,99 

40  IC-MC+l 

PXI ICI=-Tl*( X2-XI I ) i / I Z2-Z ( I,J+l) )  ♦  X  ( I ) 

PYI  IC  )=-  T  1*1  Y2-YI  JM  )  )  /(  Z2-ZI  I  ,  J*  1)  I  ♦  Y  I  J+I  J 
KCHK! 4 )  =  I 

TEST  SEGMENT  5 

49  T2=  Z I  I  +  l.J+l  )-Z0 
C=Tl*T2 
IFID)5G, 53, 59 

50  I C= I C ♦ I 

PXI  IC)=-T1*<  XI I+ll-XI I ) ) / 1 Z U+l ,J*1)-Z(I , J*1 )  H-XI I) 

PYI  1C)=Y(J<  1  ) 

KCHK I  51  =  1 
GC1 0  59 

53  II-IT2.NF.0.J)  GOTO  59 
IC=  1C  +  1 
PXI  I C  1  =  X  I  1  +  1  ) 

PYI  I  C  )  =  Y  I  J  ♦  1  ) 

KCHK I  5)  =  1 

TEST  SEGMENT  6 

59  D=T2*T3 
IFI0I60, 69,65 

60  I C= I C ♦  1 

px(ic)=-r2*(x2-xii  +  i))/iz2-zii*i,j*i))^xii  +  n 

PYI  IC)=-T2*I  Y2-YI  J+l)  )/!Z2-Z(  1*1,  J  +  m+Y  I  J*l) 

KCHK ( 6  )  =  1 
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TEST  SEGMENT  7 

69  T  1  =  T2 

T2*Z  (  I  ♦  1 . J I - Z  D 

0=Tl*T2 

IF( 0) 70, 79, 79 

70  I C= I C ♦  1 

PX<  IC)  =  X(  1*1  ) 

PYI  IC>«-T1*IY(J)-Y(  J  +  J))/IZ(I  +  1  ,J)-Zn*l,J+l))+Y<  J*l) 

KCHK ( 7 )  —  1 

TEST  SEGMENT  8 

79  D=T2*T3 

IF( D )Q0, 89, 89 

80  IC=  IC ♦  1 

PX(  I C  I  =  -  T  3*  (  X  (  I+l)-X2)/(Z(I*-l,J)-Z2)-*-X2 
PY(IC)=-T3*IY(J)-Y2)/IZ(I*l,J)-Z2»*Y2 
K CHK  ( 3 ) -  1 

89  I F (  IC.GE.2)  GOTO  90 
GOTO  201 

THIS  SFCTI0.N  OF  COOING  ORDERS  THE  DATA  SO  THAT  NO  OVERLAPPING  OR 
BACKTRACKING  OCCURS  DURING  PLOT! ING. 

90  I  F (  IC.GE  .6)  GOTO  300 
1FIKCHKI dl.NE.lI  GOTC  2CC 
DO  111  1=1,7 

IFIKCHKIL J.NE.I)  GOTO  101 

IC=  IC+1 

PXI IC )=P<( 1 ) 

PY (  I C  )  =P  Y (  1  ) 

I  C=  I C- 1 
DO  102  M=l, 1C 
P  X  (  M  )  =  P  X  (  M  ♦  1  ) 

102  PY(M)=PY(M*11 

I F ( (MCO(L,2) .EO. l).ANO.KCHK(L).EQ.l)  GOTO  200 
101  CONTINUE 

PLOT  DATA  FOR  THE  USUAL  CELL. 

2C0  CALL  PL  T  CCD  ( 1 , 0, PX II) ,P Y (1) , I C  ,C ) 

GOTO  201 

THIS  SECTION  CF  CODING  DOES  THE  ORDERING  AND  PLOTTING  IF  2  CONTOUR 

lines  pass  through  the  cell 

300  IF(Z(  I,J).GT.ZC.AND.ZII+1,J+1).GT.0.0)  GOTO  301 
1FIZ2.GT.Z0)  GOTO  3C3 

302  N  = I C ♦ 1 

F  X ( N ) =P  X ( 1 ) 

PY(N)=PY(  1) 

CALL  PL  T CCD  ( 1 , 0,  PX ( 5 ) , P Y I  5 ) , 3, C ) 

CALL  PL  T CCD  ( 1 , 0 , PX ( 2 ) ,P Y ( 2 ) , 3, 0 ) 

GOTO  311 

303  CALL  PLTCCO  ( 1, 0,PX ( 1 ) ,PY ( 1 > , 3, 0) 

CALL  PLTCCO  I l , C , PX ( A ) ,P Y ( 4  I , 3, C ) 

GOTO  31  } 

301  IFIZ2.Gr.ZG)  GOTO  3C2 
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GOTO  30  3 
310  CON TnU  = 

CLEAR  WORKING  ARRAYS  AND  INITIALIZE  FLAGS 

2C1  1C=C 

DC  *  MN=l,B 
PX( MM  »  =  : . 

pyi mni»:  . 

8  KCHK(*N>=C 
100  CONT  1  iur. 

1  CCNTINU- 

2  CCNTl'IUc 
RETURN 
END 


CONT  R  18  1 
C0NTR182 
CGNTR183 
CCNTP184 
CONT R  13 5 
CCNTR186 
C0NTR187 
C0NTR188 
C0NTR189 
CONTR  190 
CONTR 111 
CON TR  19  2 
CONTR 193 
CONTR  194 
CONT  R  195 
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